First-order phase transitions: A study through the parallel tempering method 
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We study the applicability of the parallel tempering method (PT) in the investigation of first- 
order phase transitions. In this method, replicas of the same system are simulated simultaneously 
at different temperatures and the configurations of two randomly chosen replicas can occasionally be 
interchanged. We apply the PT for the Blume-Emery-Griffiths (BEG) model, which displays strong 
first-order transitions at low temperatures. A precise estimate of coexistence lines is obtained, 
revealing that the PT may be a successful tool for the characterization of discontinuous transitions. 
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I. INTRODUCTION 

Due to the absence of exact solutions on most sys- 
tems, Monte Carlo methods play an important role not 
only in statistical physics and critical phenomena but 
also in other areas. Usually, the Metropolis [H and the 
Glauber algorithms are used to lead the system to 
the Gibbs distribution. Despite their simplicity and gen- 
erality, difhculties appears in studying the emergence of 
phase transitions when they are used to generate the mi- 
croscopic configurations. Several techniques have been 
proposed to circumvent these difficulties, such as the mul- 
ticanonical technique cluster algorithms (that work 
properly not only for reducing critical slowing down 
but also for eliminating metastability in first-order tran- 
sitions [liH J3]), the Wang-Landau method Q, simulated 
tempering 9l| , and replicas exchanges also named parallel 
tempering methods (PT) J^O, 11]. 

Special attention has been devoted to this latter ap- 
proach, due to its relative simplicity in comparison with 
other approaches and its enormous applicability for sev- 
eral systems in the framework of both statistical mechan- 
ics [12, [13, H, ^iid molecular dynamics [Tsl, fl6| . Es- 
sentially, the PT consists of simulating simultaneously a 
given set of replicas of the same system at different tem- 
peratures and, occasionally, interchanging the configura- 
tion of two randomly chosen elements of those replicas. 
This exchange between pairs of replicas allows for the 
implementation of an ergodic walk in the configuration 
space when the elements of the pair are separated by 
large free energy barriers. 

Although the PT has been widely used in several con- 
texts, an open question concerns its applicability for the 
investigation of first-order phase transitions [lo| . In fact, 
since in discontinuous transitions a gap in the energy 
might lead to a small probability of accepting exchanges 
between replicas this appears not to be a favorable sce- 
nario for PT. 

In this paper, we give a further step in this direc- 
tion by applying the PT to study and characterize first- 



order transitions. We will consider the well known spin-1 
Blume-Emery-Griffiths (BEG) model [l3|, which possess 
a rather rich phase diagram with different structures, in- 
cluding first-order transitions in the regime of low tem- 
peratures. As we shall see, the PT can be applied be- 
cause thermodynamic properties are actually described 
by continuous functions when finite systems are simu- 
lated. In fact, the discontinuity of thermodynamic prop- 
erties occur only in the thermodynamic limit. However, 
smooth curves are obtained only when one uses a dy- 
namics yielding a correct sampling of the configuration 
space 0, 0, Hi- In particular, the use of the PT allows 
for applying a new finite size procedure for the study of 
first-order phase transitions, as proposed in Ref. 0|- It is 
worth mentioning that a PT-based analysis of first-order 
transitions has recently been proposed by Neuhaus et al 
[is] . Such approach is, however, rather different from the 
one adopted here. 

This paper is organized as follows: In Sec. II we 
present the model, in Sec. Ill we describe the PT, in 
Sec. IV we discuss the numerical results, and in Sec. V 
the conclusions. 



II. MODEL 

The spin-1 BEG model is described by the following 
Hamiltonian: 



7i = — J y ^ (TiCT j 



(1) 



where ai denotes the spin variable of the z-th site of the 
lattice which assumes the values ±1 and and the sums 
run over the nearest neighbor spins on a d— dimensional 
lattice with V = L"^ sites. Parameters J, K are the 
nearest-neighbor spin couplings and D is the quadrupole 
moment. We have two order parameters defined as fol- 

lows: q = {j:ti'^f)/y and m = (ELi '^»)/^- The 
BEG model will be consider for a square lattice and pe- 
riodic boundary conditions. 
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III. PARALLEL TEMPERING METHOD 

In the parallel tempering method (PT), configurations 
at high temperatures are used to perform an ergodic walk 
in low temperatures. To this end, we simulate, for fixed 
values of D, a set of N replicas in the interval of temper- 
atures {Ti, ...,Tn}, where Ti and Tjv are extreme tem- 
peratures. 

The dynamics is composed of two parts. In the first 
part, each one of the N replicas are simulated accord- 
ing to the Metropolis algorithm. For the i~th replica 
a given site k of the system is chosen at random and 
we select, with equal probability, one of the two other 
possible spin values a'^ such that cr^ ^ ak- The spin 
variable ak is then replaced with according to the 
Metropolis prescription: pk — min{l, exp(— /SAT^)} 
where AH = n{a'^) - n{ak) and /3 = l/ksT. In the 
second part of the dynamics, the PT is implemented. 
After a given number of Monte Carlo steps, the ex- 
change of configurations of two replicas at the tem- 
peratures Ti and Tj are performed with the probabil- 
ity pij = min{l,exp[(/3i - j3j){TL{ai) - 7i((Tj)]}, where 
Tj > Ti, j = i + S, and 5 denotes the "distance" between 
two arbitrary replicas. The probability pij depends on 
{(3i — (3j ) and for this reason the performance of method 
will depend on the "distance" between the replicas. If 
the difference is large enough exchanges will be hardly 
performed and the PT will not provide any improvement 
in the results. 

In this paper, we adopt two independent procedures 
to choose the interval of temperatures. In the first one, 
the distance between adjacent temperatures obey a ge- 
ometric progression. Some authors have shown jlOj, |20| 
that while this procedure works well when specific heat of 
the system is about constant, at the emergence of a phase 
transition, when the specific heat diverges, its efficiency is 
reduced. For this reason, we adopted a second procedure, 
that consists in distributing temperatures in regular in- 
tervals between Ti and T^ for a given small size system. 
By increasing L, we introduce additional temperatures 
between T and T^+i. This procedure is necessary be- 
cause the exchange probability in general decreases as L 
increases. We have verified that both procedures lead to 
the same results, within of the statistical errors. 

Concerning the replicas exchanges we also consider 
exchanges between nonadjacent sites. This is imple- 
mented in this work by allowing 5 to range in the in- 
terval 5 = 1,..,6. As it will be shown, although non- 
adjacent exchanges have been less studied EBli be- 
cause the probability of performing a given exchange de- 
creases when 5 increases, they have revealed to be essen- 
tial mechanisms in eliminating hysteretic effects. 

IV. NUMERICAL RESULTS 

We have simulated three different values of K/J, 
given by K/J =0, 3, and 3.3. Note that the first 




FIG. 1: Order parameter g as a function of D for K/J = 3, 
T — 1.5 and L = 30 obtained from parallel tempering (sym- 
bol x) and cluster algorithms (circles). Squares correspond 
to data obtained from parallel tempering with exchanges only 
between nearest-neighbor replicas. In the inset, circles and 
triangles refer to the Metropolis algorithm, whereas the sym- 
bol X refers to the parallel tempering. 

case {K/J — 0) corresponds to the well known Blume- 
Capel model. Replicas are distributed in the intervals 
Ti = 1.5 < T < 2.2 = Tn, for K/J = 3 and 3.3, and 
Ti = 0.4 < T < 0.62 Tn, for K/J = 0. We have 
simulated system with size L ranging from L = 10 up 
to L = 40 and we considered 8 x 10^ Monte Carlo steps 
to evaluate the appropriate quantities after equilibrating 
the system. For all values of K/J considered here, the 
system displays two ferromagnetic phases (rich at spins 
-I- and — ) for small values of D. Also, a paramagnetic 
phase (rich at spins 0) takes place for high values of D. A 
strong first-order phase transition between the ferromag- 
netic and paramagnetic phase occurs for a given value of 
D that depends on K/J and T. 

The first inspection about the applicability of the PT 
for first-order transitions is shown in the inset of Fig. [U 
where we compare the PT results with those obtained by 
using only the Metropolis algorithm. By simulating only 
with the Metropolis algorithm, the system gets trapped 
in metastable states and even after 8 x 10^ MC steps it 
does not undergo a transition to the stable phase. This 
effect does not occur when we use the PT with nonlo- 
cal exchanges, since the system becomes able to pass 
from one phase to the other. The efficiency of the PT 
is also corroborated by the agreement with results ob- 
tained from cluster algorithms Q , where a smooth curve 
is obtained for the order parameter. However, as it was 
mentioned previously, when one considers only exchanges 
of configurations between nearest-neighbor replicas, hys- 
teresis are still present, as showed in FiglTJ 

The role of non-local exchanges is analyzed in more 
details by considering the time evolution of thermody- 
namic properties at the phase coexistence. In Fig. [2] 
we plot, for a single run, the order parameter q start- 
ing from two different initial configurations for K/J = 3, 
T = 1.5 and L = 20. In the inset of each graph, we 
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FIG. 2: Time evolution of the order parameter q for a sin- 
gle run starting from two independent initial configurations 
simulated with (a) the Metropolis algorithm and (6) the PT, 
for L = 20, T = 1.5, D = 8.0, and K/J = 3. In the insets 
the time evolution of the total energy per volume u is given 
for those initial configurations. In contrast with the PT, until 
M = 6 X 10'* MC steps the Metropolis algorithm provides a 
nonergodic simulation. 
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FIG. 3: In graph (a) we show the time evolution of the or- 
der parameter q simulated by cluster algorithm and PT with 
exchanges between i and its i -\- 5th next neighbor replica 
{5 = 1,2,3 and 6) for L = 20 and NR = 500 independent 
runs. In graph (b) it is shown the mean probability p* versus 
T for different 5. 



plot the time evolution of the total energy per volume 
u for the same initial configurations. In contrast with 
the PT, until M = 6 x 10'' MC steps, the simulation 
is not ergodic when the system is simulated with the 
Metropolis algorithm. Next, in Fig. [3] (a) the time evo- 
lution of the system simulated via PT with local and 
non-local exchanges is comapared with the results pro- 
vided by cluster algorithms. Note that for 5 > 2 and 
M > 3 X 10^ MC steps, the time evolution the PT sim- 
ulation for q converges to (7 « 2/3 (as will be explained 
later), in agreement with cluster algorithm simulations. 
A similar behavior is obtained in all cases for the quantity 
m. In Fig. [3] (6) shows the exchange mean probability 
p* = (min{l,exp[(A - /3j)(H(a,) - H(a,)]}) [H as a 
function of T for different distances 5 between replicas 
and L = 20. Except for 5 — \, the minimum in p* oc- 
curs at T « 1.95, indicating the coexistence between the 
ferromagnetic phases, paramagnetic rich at spins and 
a disordered phase, that takes place in the limit of high 
temperatures |T^ . Our results show that, although non- 
local exchanges are performed less frequently than local 
ones, they are fundamental for ensuring an ergodic simu- 
lation of the system. Next, we will describe the method- 
ology employed in determining coexistence lines. Their 
location will be derived from finite size analysis for both 
the order parameter q and the susceptibility xt- 

Although a discontinuous phase transition is character- 
ized by a jump in the order parameter, the discontinuity 
takes place only in the thermodynamic limit. For finite 
systems, not only the order parameter, but also other 
quantities are described by continuous functions 0, [I]. 
In such case, the behavior of pWsical quantities scales 
with the volume of the system [51, [2^. In Fig. 01 the 
order parameter q is shown as a function of D for several 



values of L. 

Although isotherms present strong dependence on the 
system size, they intersect one another at the point 
D = = 8.0000(1) and g « 2/3. As it was explained 
m Refs. 0, [11 by' means of two different reasonings, 
the point where all isotherms cross is independent of the 
lattice size. This can be understood recalling that in 
the regime of low temperatures two ferromagnetic phases 
(g « 1) coexist with a paramagnetic phase rich at spins 
zero {q K, Q) dX D — Dq, yielding q w 2/3 for all sys- 
tem sizes. Therefore, the crossing point can be used as 
a criterion to estimate the transition point. As it will 
be shown later, the estimate of Dg agrees very well with 
the value obtained from finite size analysis for the 
susceptibility xt- In Fig. IHb), we describe the collapse 
of all data by the expression y* = {D — D'o)L^ confirm- 
ing the dependence on the volume. At low temperatures, 
the relation between q with the system size L and D is 
expressed by the following equation 0, [2^ 



q = 



b + ce"°^ 



1 + de- 



(2) 



where a, 6, c and d are fitting parameters and z = D~Dq. 
In Fig. [4] (a), continuous lines correspond to the fittings 
proposed by Eq. The parameter d scales with the 
volume, as shown in Fig. (He). In the thermodynamic 
limit L — > 00, while the quantity d diverges the order 
parameter q does not. According to Eq. ([2]), in the ferro- 
magnetic phase, which occurs in the region D — Dq < 0, 
we have that q — > c/d as L — > 00. On the other hand, in 
the paramagnetic phase, which appears for D — Dq > 0, 
g — > 6 as i — > 00. For D = Dq, we have a jump in q, 
signing a discontinuous phase transition. 

In the second analysis, we determine the transition 
point by examinating the susceptibility xt — PLp'{{<f') — 
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FIG. 4: Order parameter per volume q versus D for several 
values of the system size L for K/ J = 3 and T = 1.5. Con- 
tinuous lines correspond to the fittings defined by Eq. ((Sjl. 
In (6) we have a collapse of all data by using the relation 
y* — [D — Dq)L^ . In (c) we have the log-log plot for the 
quantity a as a function of L. The straight line has slope 
2.00(1). 



(q)^). Increasing D towards the coexistence line, one ob- 
serves a sharp peak in xt at D*^ for all system sizes, 
as shown in Fig. El^a). The deviation between D]^ and 
its asymptotic value D*^ decays as in a first-order 
transition [l^, . Our results satisfy this asymptotic re- 
lation, as it can be seen in Fig. EK^). From this law, we 
have obtained the extrapolated value Z?^ — 8.0000(1), 
which agrees with the estimate Dq obtained previously 
and also agrees with the result D = 8.0000(1), obtained 
from a cluster algorithm for the BEG model [7]. In Fig. [5] 
(c) we observe that all curves coalesce to x* = Xt/ L'^ and 
y* = {D1 — D'^)L^, confirming once again the scaling 
with the volume. 

It is worth emphasizing that when one uses only the 
Metropolis algorithm to generate the configurations, nei- 
ther the crossing among isotherms nor accurate finite size 
analysis for smooth curves become possible, due to the 
presence of hysteresis effects, as it can be seen in Fig. [T] 

In Figs. and [71 we repeat, for K/J = 3.3 and 
T ~ 1.5, both analysis presented above for determin- 
ing phase coexistence. From the first procedure, where 
all isotherms are to be fitted by Eq. the crossing is 
given by g « 2/3 and = 8.6032(1). This estimate 
agrees with the value = 8.6033(1) obtained from fi- 
nite size analysis for the quantity xt, as showed in Fig. 
[71 These estimates, both obtained by using the PT, are 
in good accordance with the value D = 8.6032, obtained 
by Rachadi and Benyoussef, from cluster algorithms 6]. 



In the last analysis, we show in Fig.[51numerical results 
for K/J — considering T — 0.4. Fitting all isotherms 
with Eq. ^ , the intersection point turns out to be given 
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FIG. 5: Susceptibility xt versus D for several values of the 
system size L, K/J = 3, and T — 1.5. In (6), we plot the 
value of D = Z)|, in which xt is maximum, as a function of 
. In (c) we have a collapse of all data using the relations 

X* = xt/L^ and y* = (D - D'^)L^. 
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FIG. 6: Order parameter per volume q versus D for several 
values of the system size L for K/J = 3.3 and T = 1.5. 
Continuous lines stand for the fittings defined by Eq. 
In (6) we have a collapse of all data using the relation y* = 
{D — Dq)L^ . In (c) we have the log-log plot of a versus L. 
The straight line has slope 2.00(1). 



by g « 2/3 and = 1.9968(1). The collapse of data 
using this estimate of Dq confirms again the adequacy 
of this procedure for the estimation of transition point. 
Repeating this procedure for T = 0.5, we verify that all 
isotherms cross the abscissa Dq = 1.9879(1), which is 



in fair agreement with the estimates T = 0.499(3) 
D = 1.992, predicted by Wang-Landau method [2l| 



and 
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FIG. 7: Susceptibility xt versus D for several values of the 
system size L, K/ J = 3.3, and T = 1.5. In (6), we plot the 
value of D = Z)|, in which xt is maximum, as a function of 
. In (c) we have a collapse of all data using the relations 
X* = Xr/i' and y' = {D - 0*^)1". 
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CONCLUSION 



In this paper, we have apphed the parallel tempering 
-method (PT) for the study of first-order transitions. We 
nave considered different regions of the phase diagram 
0.004 _20. 008 of the BEG model, for which usual Metropolis algorithm 
leads to strong hysteresis at the phase coexistence, pro- 
j^iding no reliable estimates of the coexistence lines. On 
the other hand, by using the PT it was possible to cir- 
cumvent the free energy barriers and as a consequence 
Hiysteretic effects were eliminated. All results obtained 
via PT allowed us to locate the transition points pre- 
cisely by means of two techniques, whose estimates agree 
j^gith those obtained from other procedures, such as clus- 
ter algorithms and, in one case, with the Wang-Landau 
method. Although the agreement between results ob- 
tained from PT and cluster algorithms have been shown 
to be very well, cluster algorithms are more specialized, 
since each model requires a specific cluster algorithm that 
takes into account the appropriate transitions. On the 
other hand, PT is general and can be used, in principle, 
for any system. We remark that more studies of first- 
order transitions using the parallel tempering are still 
required, in order to have a more comprehension of its 
performance. 
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FIG. 8: Order parameter per volume q versus D for several 
values of system size L, for K/ J = and T = 0.4. Continuous 
lines stand for the fittings defined by Eq. ([2]). In the inset, we 
have a collapse of all data using the relation y* — {D — Dg )L^ . 
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